Gradient Forest - Genomic offset predictions

Introduction

R code based on the github repository associated with Fitzpatrick et al. (2021).

Data

Genomic data

Code

# Population-based allele frequencies
# ===================================
geno <- read.csv(here("data/DryadRepo/ImputedGenomicData_AlleleFrequencies_withoutmaf.csv"),
                     row.names = 1)

# SNP sets
# ========
snp_sets <- readRDS(here("outputs/list_snp_sets.rds"))

Climatic data

Code

# Set of climatic variables
# =========================
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))


# Climatic data
# =============
# we scale the past and future climatic data at the location of the populations
# with the parameters (mean and variance) of the past climatic data.
source(here("scripts/functions/generate_scaled_clim_datasets.R"))
clim_dfs <- generate_scaled_clim_datasets(clim_var,clim_ref_adj = FALSE)

Running GF models

Code

snp_sets <- lapply(snp_sets, function(snp_set) {


# Warning! Important to sort the SNP names, otherwise
  # the colors in the species (allele) cumulative plots 
  # do not correspond to the right alleles
geno_sub <- geno %>% dplyr::select(all_of(sort(snp_set$set_snps)))

snp_set$gf_mod <- gradientForest(data.frame(clim_dfs$clim_ref[,-1], geno_sub), 
                                 predictor.vars=clim_var, 
                                 response.vars=colnames(geno_sub), 
                                 corr.threshold=0.5, 
                                 ntree=500, 
                                 trace=T)

return(snp_set)

})

Warning in randomForest.default(x = X, y = spec_vec, maxLevel = maxLevel, : The response has five or fewer unique values.  Are you sure you want to do regression?

Number of alleles for which the climatic variables have some predictive power:

Code

lapply(snp_sets, function(x){
  tibble("SNP set"=x$set_name,
             "Nb of alleles with predictive power"=x$gf_mod$species.pos.rsq)
}) %>% 
  bind_rows() %>% 
  kable_mydf()

SNP set	Nb of alleles with predictive power
Candidate SNPs	69
Candidate SNPs (with correction)	45
Control SNPs	46

We generate some plots to evaluate the GF models (stored in the files GFplots_[SNP set code].pdf):

• Predictor overall importance plots. They show the mean accuracy importance and the mean importance weighted by SNPs \(\mathcal{R}^2\).

• Splits density plots. They show the binned split importance and location on each gradient (spikes), kernel density of splits (blacklines), of observations (red lines) and of splits standardised by observations density (bluelines). Each distribution integrates to predictor importance. These show where important changes in the abundance of multiple alleles are occurring along the gradient; they indicate a composition change rate.

• Species (in our case alleles) cumulative plots. They show, for each SNPs, the cumulative importance distributions of splits improvement scaled by \(\mathcal{R}^2\) weighted importance, and standardised by density of observations. These show cumulative change in the presence of individual allele, where changes occur on the gradient, and the alleles changing most on each gradient.

• Predictor cumulative plots. They show, for each predictor, the cumulative importance distributions of splits improvement scaled by \(\mathcal{R}^2\) weighted importance, and standardised by density of observations, averaged over all SNPs. These show cumulative change in overall allelic composition, and where changes occur on the gradient.

• \(\mathcal{R}^2\) measure of the fit of the random forest model for each SNPs.

The code to generate those plots comes from ‘Example analysis of biodiversity survey data with R package gradientForest’ by C. Roland Pitcher, Nick Ellis and Stephen J. Smith (pdf available here).

Code

# ==============================
# Functions to make the GF plots
# ==============================

# splits density plot
make_split_density_plot <- function(x){

  plot(x,plot.type="S",
       imp.vars=names(importance(x)),
       leg.posn="topright",
       cex.legend=1,
       cex.axis=1, 
       cex.lab=1.2,
       line.ylab=0.9,
       par.args=list(mgp=c(1.5, 0.5,0),mar=c(3.1,1.5,0.1,1),omi = c(0.1, 0.3, 0.1, 0.1)))
  
  }

# predictor cumulative plot
make_predictor_cumulative_plot <- function(x){
  
  plot(x, plot.type="C",
       imp.vars=names(importance(x)),
       show.species=F,common.scale=T,
       cex.axis=1, 
       cex.lab=1.5,
       line.ylab=1,
       par.args=list(mgp=c(1.5, 0.5,0),mar=c(2.5,1,0.1,0.5),omi=c(0, +0.3,0,0)))}


# species cumulative plot
make_species_cumulative_plot <- function(x){
  
  plot(x,
     plot.type="C",imp.vars=names(importance(x)), 
     show.overall=F,legend=T,leg.posn="topleft", 
     leg.nspecies=10,
     cex.lab=1,
     cex.legend=1, 
     cex.axis=1,
     line.ylab=1,
     par.args=list(mgp=c(1.5, 0.5,0),mar=c(2.5,1,0.1,0.5),omi=c(0, +0.3,0,0)))
  }


# R2 measure of the fit of the random forest model for each species
make_performance_plot <- function(x, horizontal=FALSE){

    old.mar<-par()$mar
    par(mfrow=c(1,1),mar=old.mar+c(0,0,0,0))
    Ylab <- expression(R^2)

    perf <- importance(x, type="Species")
    n <- length(perf)
    if (horizontal)
      plot(perf, 1:n, las = 2, pch=19, axes=F, xlab="", ylab="")
    else
      plot(1:n, perf, las = 2, pch=19, axes=F, xlab="", ylab="")
    axis(labels= names(perf),side=1+horizontal,at=1:n, cex.axis=0.7, padj=0,las=2)
    axis(side=2-horizontal, cex.axis=1)
    mtext(Ylab,side=2-horizontal,line=2)
    title("Overall performance of random forests over loci")
    abline(h = 0, lty = 2)
    box()
    par(mar=old.mar)

}

# ===============================================
# Generate the GF plots for the Github repository
# ===============================================
lapply(snp_sets, function(snp_set){

pdf(here(paste0("figs/GF/GFplots_",snp_set$set_code,".pdf")), width=12,height=8)

# Overall importance plot
plot(snp_set$gf_mod, plot.type="Overall.Importance")

# splits density plot
make_split_density_plot(x=snp_set$gf_mod)

# species cumulative plot
make_species_cumulative_plot(x=snp_set$gf_mod)

# predictor cumulative plot
make_predictor_cumulative_plot(x=snp_set$gf_mod)

# R2 measure of the fit of the random forest model for each species
make_performance_plot(x=snp_set$gf_mod)

dev.off()

})

# =======================================================
# Generate the GF plots for the Supplementary Information
# =======================================================
lapply(snp_sets[c(1,3)], function(snp_set){

# Overall importance plot
pdf(here(paste0("figs/GF/GFplots_OverallImportance_",snp_set$set_code,"_SI.pdf")), width=8,height=5)
plot(snp_set$gf_mod, plot.type="Overall.Importance")
dev.off()

# Splits density plot
pdf(here(paste0("figs/GF/GFplots_SplitDensityPlot_",snp_set$set_code,"_SI.pdf")), width=8,height=5)
make_split_density_plot(x=snp_set$gf_mod)
dev.off()

# Allele cumulative plot
pdf(here(paste0("figs/GF/GFplots_AlleleCumulativePlot_",snp_set$set_code,"_SI.pdf")), width=7,height=7)
make_species_cumulative_plot(x=snp_set$gf_mod)
dev.off()

# predictor cumulative plot
pdf(here(paste0("figs/GF/GFplots_PredictorCumulativePlot_",snp_set$set_code,"_SI.pdf")), width=7,height=7)
make_predictor_cumulative_plot(x=snp_set$gf_mod)
dev.off()

# R2 measure of the fit of the random forest model for each species
p <- tibble(snp=names(sort(snp_set$gf_mod$result)),importance=sort(snp_set$gf_mod$result)) %>% 
  ggplot() +
  geom_point(aes(y=reorder(snp, importance),x=importance)) +
  ylab("") +
  xlab(expression(R^2)) +
  theme_bw() 

ggsave(p, filename = here(paste0("figs/GF/GFPlots_AlleleImportance_",snp_set$set_code,"_SI.pdf")), width=4, height=10)


})

GO predictions

Code

snp_sets <- lapply(snp_sets, function(snp_set){
  
snp_set$go <- lapply(clim_dfs$clim_pred, function(clim_pred){

ref_pred <- predict(snp_set$gf_mod) # predictions under current climates
fut_pred <- predict(snp_set$gf_mod, as.data.frame(clim_pred[,clim_var])) # predictions under future climates


lapply(1:nrow(ref_pred), function(x, ref_pred, fut_pred){
    as.numeric(pdist(ref_pred[x,],  fut_pred[x,])@dist)}, fut_pred=fut_pred, ref_pred=ref_pred) %>% 
  unlist()

})
return(snp_set)
})

Relationship with Euclidean distance

Code

source(here("scripts/functions/make_eucli_plot.R"))

# Calculate the Euclidean climatic distance
list_dist_env <- clim_dfs$clim_pred %>% lapply(function(clim_pred){
  
Delta = clim_dfs$clim_ref %>% dplyr::select(any_of(clim_var)) - clim_pred %>% dplyr::select(any_of(clim_var)) 
dist_env = sqrt( rowSums(Delta^2) )

})

# Main gene pools (for the figures)
gps <- readRDS(here("data/GenomicData/MainGenePoolPopulations.rds")) %>%  arrange(pop)

Code

par(mfrow=c(1,2))

lapply(snp_sets, function(x) {

lapply(names(list_dist_env), function(gcm){
  
make_eucli_plot(
  X = list_dist_env[[gcm]],
  Y = x$go[[gcm]],
  colors = gps$color_main_gp_pop,
  color_names = gps$main_gp_pop,
  ylab = "GF genomic offset",
  legend_position="topleft",
  plot_title = paste0(x$set_name," - ", gcm))

})
}) 

Code

# We generate scatter plots for the Supplementary Information.
# ============================================================

# Axis limits
# ===========
max_go <- lapply(snp_sets[c(1,3)], function(z){
  z$go %>% unlist()
}) %>% unlist() %>% max()

range_eucli <- list_dist_env %>% unlist() %>% range()

# Run the function
# ================
lapply(snp_sets[c(1,3)], function(set_i) {

p <- lapply(names(list_dist_env), function(gcm){
  
make_ggscatterplot(
  x = list_dist_env[[gcm]],
  y = set_i$go[[gcm]],
  title=gcm,
  range_eucli = range_eucli,
  max_go = max_go)

})

# remove y-labels to graphs in the second column
p[[2]] <- p[[2]] + ylab("")
p[[4]] <- p[[4]] + ylab("")


# remove x-labels to graphs in the second and third rows
p[[1]] <- p[[1]] + xlab("")
p[[2]] <- p[[2]] + xlab("")
p[[3]] <- p[[3]] + xlab("")

p[[6]] <- get_legend(p[[1]])

for(i in 1:5){p[[i]] <- p[[i]]  +  theme(legend.position = "none")} 


plot_grid(plotlist=p, nrow = 3) %>% 
  ggsave(here(paste0("figs/GF/ScatterPlotEucliDistance_",set_i$set_code,".pdf")), 
         .,
         width=7,
         height=8,
         device="pdf")

})

Comparing GO predictions

We look at the correlation across the different genomic offset predictions at the location of the populations, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code

par(mfrow=c(1,3))

lapply(names(snp_sets[[1]]$go), function(gcm){
  
lapply(snp_sets, function(x) x$go[[gcm]]) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', 
           diag = FALSE,mar=c(0,0,2,0),
           title=gcm,
           number.cex=2,tl.cex=1.5)
  
})

Maps

With spatial points

Code

# Function to make the genomic offset maps
source(here("scripts/functions/make_go_map.R"))

# Population coordinates
pop_coord <-  readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_noADJ.rds")))[[1]]$ref_means %>% dplyr::select(pop,longitude,latitude)

# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(x) {
lapply(names(list_dist_env), function(gcm){
x$go[[gcm]]
}) %>%  unlist()
}) %>%  unlist() %>% range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0


# Generate the maps for each set of SNPs and each GCM
lapply(snp_sets, function(x) {

go_maps <- lapply(names(list_dist_env), function(gcm){
  
make_go_map(dfcoord=pop_coord,
            snp_set = x,
            gcm=gcm,
            ggtitle=gcm,
            go_limits = go_limits,
            point_size = 3)

})

legend_maps  <- get_legend(go_maps[[1]])

go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))

go_maps$legend_maps <- legend_maps

go_maps <-plot_grid(plotlist=go_maps)

# =====================================================
# We save the figures for the Supplementary Information
# =====================================================
if(x$set_code=="cand"){ 
  ggsave(here("figs/GF/GOMaps_PopLocations_CandidateSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")
  
  } else if(x$set_code=="control"){
    ggsave(here("figs/GF/GOMaps_PopLocations_ControlSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")
    }


# =========
# Add title
# =========
title <- ggdraw() + 
  draw_label(
    x$set_name,
    fontface = 'bold',
    x = 0,
    hjust = 0
  ) +
  theme(plot.margin = margin(0, 0, 0, 7))

# merge title and plots
plot_grid(
  title, go_maps,
  ncol = 1,
  rel_heights = c(0.1, 1))
})

For each GCM, we attribute the value 1 to the top five populations with the highest genomic offset and we attribute the value 0 to the other populations. We then count the number of 1 for each population, which gives the table and map below:

Code

source(here("scripts/functions/make_high_go_pop_maps.R"))

high_go_pops <- make_high_go_pop_maps(pop_coord=pop_coord,
                                      list_go = snp_sets$cand$go,
                                      ggtitle="GF",
                                      nb_id_pop = 5) # number of selected populations

saveRDS(high_go_pops, file = here("outputs/GF/high_go_pops.rds"))

high_go_pops[[1]] %>% kable_mydf

pop	longitude	latitude	GFDL-ESM4	IPSL-CM6A-LR	MPI-ESM1-2-HR	MRI-ESM2-0	UKESM1-0-LL	sum_go
PUE	-6.63	43.55	0	0	0	0	0	0
LEI	-8.96	39.78	0	0	0	0	0	0
PIA	9.46	42.02	0	0	0	0	0	0
QUA	-0.36	38.97	0	0	0	0	0	0
OLB	-0.62	40.17	0	0	0	0	0	0
LAM	-6.22	43.56	0	0	0	0	0	0
COC	-4.50	41.25	0	0	0	0	0	0
CUE	-4.48	41.37	0	0	0	0	0	0
BON	-1.66	39.99	0	0	0	0	0	0
TAM	-5.02	33.60	0	0	0	0	0	0
CAR	-4.28	41.17	0	0	0	0	0	0
MAD	-5.23	35.18	0	0	0	0	0	0
SAL	-3.06	41.83	0	0	0	0	0	0
CAS	-6.98	43.50	0	0	0	0	0	0
CAD	-6.42	43.54	0	0	0	0	0	0
BAY	-2.88	41.52	0	0	0	0	0	0
PLE	-2.34	47.78	0	0	0	0	0	0
HOU	-1.15	45.18	0	0	0	0	0	0
STJ	-2.03	46.76	0	0	0	0	0	0
SIE	-6.49	43.53	0	0	0	0	0	0
ORI	-2.35	37.53	0	0	0	0	0	0
SAC	-8.36	42.12	0	0	0	0	0	0
PIE	9.04	41.97	0	0	0	0	0	0
VER	-1.09	45.55	0	0	0	0	0	0
ARN	-5.12	40.19	0	0	0	0	0	0
CEN	-4.49	40.28	0	0	1	0	0	1
VAL	-4.31	40.52	0	0	1	0	0	1
OLO	-1.83	46.57	0	0	0	1	1	2
MIM	-1.30	44.13	1	1	0	0	0	2
PET	-1.30	44.06	1	1	0	0	0	2
COM	-3.95	36.83	0	1	0	1	1	3
SEG	-8.45	42.82	1	0	1	1	1	4
ALT	-6.49	43.28	1	1	1	1	1	5
ARM	-6.46	43.30	1	1	1	1	1	5

Code

high_go_pops[[2]]

With rasters

We project the genomic offset predictions for the set of candidate SNPs and the mean genomic offset across the five GCMs.

Code

# Extract scaling parameters, i.e. mean and variance  
scale_params <- lapply(clim_var, function(x){
    
    vec_var <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_noADJ.rds")))[["ref_1901_1950"]]$ref_means[,x] %>% pull()
    
    list(mean = mean(vec_var),
         sd = sd(vec_var))
    
  }) %>% setNames(clim_var)

Code

# Buffer for the maps (maritime pine distribution)
range_buffer = shapefile(here('data/Mapping/PinpinDistriEUforgen_NFIplotsBuffer10km.shp'))

# We project the genomic offset only for the candidate SNPs
snp_set <- snp_sets[["cand"]]

# We load the rasters with the climates of the reference period
path <- here("data/ClimaticData/ClimateDTRasters/1km_1901-1950_Extent-JulietteA/")
ref_rasts <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>%
  mask(range_buffer) 

# We extract the climatic values in a data frame
clim_ref_df <- ref_rasts %>% 
  rasterToPoints() %>% 
  as.data.frame()

# Scale the climatic variables with the initial scaling parameters
for(i in clim_var){clim_ref_df[,i] <- (clim_ref_df[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}

Code

# GF predictions under current climates
ref_pred <- predict(snp_set$gf_mod,clim_ref_df[,clim_var])

df <- lapply(names(clim_dfs$clim_pred), function(gcm){ # for each GCM
  
# Rasters with future climates
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
clim_fut_df <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>% 
  mask(range_buffer) %>% 
  rasterToPoints() %>% # extract the climatic values at each spatial points
  as.data.frame()

# Scale the climatic variables with the initial scaling parameters
for(i in clim_var){clim_fut_df[,i] <- (clim_fut_df[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}

# GF predictions under future climates
fut_pred <- predict(snp_set$gf_mod, clim_fut_df[,clim_var]) 

# Calculate the genomic offset
GO <- lapply(1:nrow(ref_pred), function(x, ref_pred, fut_pred){
    as.numeric(pdist(ref_pred[x,],  fut_pred[x,])@dist)}, fut_pred=fut_pred, ref_pred=ref_pred) %>% 
  unlist()

clim_ref_df %>%
  dplyr::select(x,y) %>% 
  as_tibble() %>% 
  mutate(GO = GO,
         gcm = gcm)

}) %>% 
  bind_rows %>% 
  pivot_wider(values_from = GO, names_from = gcm) %>% 
  mutate(mean_GO = rowMeans(dplyr::select(.,-c(x:y))))


df %>% saveRDS(here("outputs/GF/go_pred_rasters.rds"))

Code

# Map options
# ===========
point_size = 2
x_limits = c(-10, 15)
y_limits = c(31, 52)

# Country borders
world <- ne_countries(scale = "medium", returnclass = "sf")

# Load the mean GO projections
df <- readRDS(here("outputs/GF/go_pred_rasters.rds"))

p <- ggplot(data=df) + 
  geom_sf(data = world, fill="gray98") + 
  scale_x_continuous(limits = x_limits) +
  scale_y_continuous(limits = y_limits) + 
  geom_raster(aes(x = x, y = y, fill = mean_GO), alpha = 1) + 
  scale_fill_gradient2(low="blue", mid= "yellow", high="red",
                       midpoint=(max(df$mean_GO)-min(df$mean_GO))/2,
                       limits=c(min(df$mean_GO),max(df$mean_GO)),
                       name = "Genomic offset from GF") +
  xlab("Longitude") + ylab("Latitude") +
  theme_bw(base_size = 11) +
  theme(panel.grid = element_blank(), 
        plot.background = element_blank(), 
        panel.background = element_blank(), 
        legend.position = c(0.8,0.15),
        legend.box.background = element_rect(colour = "gray80"),
        legend.title = element_text(size=10),
        strip.text = element_text(size=11))

p %>% ggsave(here("figs/GF/GOmeanProjections_CandSNPs.pdf"),., width=8,height=8, device="pdf")

Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.

Code

p

Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
ℹ Consider using `geom_tile()` instead.

Validation - NFI plots

Predicting GO in the NFI plots

Warning

The NFI climatic datasets used for GO predictions have to be scaled with the scaling parameters used for GO estimation (i.e., the mean and variance of the climatic variables at the location of the studied populations under the reference climate).

Code

# Load the climatic data of the NFI plots.
nfi_clim <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))
  
# Keep only the climatic variables of interest and scale the climatic data
# Careful here! We have to scale the climatic variables of the reference period
# with the scaling parameters used for estimating GO!
source(here("scripts/functions/generate_scaled_nfi_clim_datasets.R"))
nfi_dfs <- generate_scaled_nfi_clim_datasets(clim_var, clim_ref = nfi_clim$clim_ref, clim_pred = nfi_clim$clim_survey)
#generate_scaled_clim_datasets(clim_var, clim_ref = nfi_clim$clim_ref, clim_pred = nfi_clim$clim_survey)

# calculate the genomic offset for the NFI plots
snp_sets <- lapply(snp_sets, function(snp_set){
  
ref_pred <- predict(snp_set$gf_mod, as.data.frame(nfi_dfs$clim_ref[,clim_var])) # predictions under reference-period climates
fut_pred <- predict(snp_set$gf_mod, as.data.frame(nfi_dfs$clim_pred[,clim_var])) # predictions under climates during survey period

snp_set$go_nfi <- lapply(1:nrow(ref_pred), function(x, ref_pred, fut_pred){
    as.numeric(pdist(ref_pred[x,],  fut_pred[x,])@dist)}, fut_pred=fut_pred, ref_pred=ref_pred) %>% 
  unlist()


return(snp_set)
})

# checking missing data
# lapply(snp_sets, function(x) sum(is.na(x$go_nfi)))

# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(snp_set) snp_set$go_nfi) %>%  unlist() %>% range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0


# map genomic offset predictions in the NFI plots 
lapply(snp_sets, function(x) {
  
  p <- make_go_map(
  dfcoord= readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))[[1]] %>% dplyr::select(contains("ude")), 
  snp_set = x,
  type="NFI",
  point_size = 0.5,
  go_limits = go_limits,
  legend_position = c(0.85,0.15),
  legend_box_background = "gray80",
  y_limits = c(35, 51))
  
  
# Figure for the SI
# =================
if(x$set_code=="cand"){ 
  p_SI <- p + theme(plot.title = element_blank())
  ggsave(here("figs/GF/NFI_GOmap_CandidateSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf")
  
  } else if(x$set_code=="control"){
    p_SI <- p + theme(plot.title = element_blank())
    ggsave(here("figs/GF/NFI_GOmap_ControlSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf")
    }

# Show maps in the Quarto document
# ================================
p
  
  })

We look at the correlation across the different genomic offset predictions in the NFI plots, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code

lapply(snp_sets, function(x) x$go_nfi) %>% 
  as_tibble() %>%
  cor() %>% 
  corrplot(method = 'number',type = 'lower', diag = FALSE,mar=c(0,0,2,0),
               number.cex=2,tl.cex=1.5)

How to explain high GO values in some of the NFI plots?

Interestingly, GO predictions for some of the NFI plots have considerably higher values than GO predictions at the location of the populations under future climates (section Section 4) or GO predictions projected across the species distribution (section Section 4.3.2) under future climates.

We can visualize these differences with the graph below that shows the distribution of GO predictions projected across the species range under future climates (i.e., with the 5 GCMs) and GO predictions at the location of the NFI plots during the survey period.

Code

readRDS(here("outputs/GF/go_pred_rasters.rds")) %>% 
  dplyr::select(-x,-y) %>%
  dplyr::select(-mean_GO) %>% # we do not show the mean GO across GCMs
  set_colnames(paste0("Species range - Future climate with ",colnames(.))) %>% 
  pivot_longer(cols = everything(), names_to = "projection") %>% 
  bind_rows(tibble(projection="NFI plots - Climate during the inventory period",value=snp_sets$cand$go_nfi)) %>% 
  
  ggplot(aes(x=value, group=projection, fill=projection)) +
  geom_density(adjust=1.5, alpha=.4) +
  theme_bw() +
  xlab("Genomic offset predictions") +
  ylab("Density") +
  theme(legend.position = c(0.7,0.8),
        legend.title = element_blank())

How can we explain these differences among the predictions?

Climatic differences

To understand, we first look at the climatic differences (1) at the location of the studied populations, (2) across the species range and (3) at the location of the NFI plots. We compare their reference climate, their climate for GO predictions (i.e., climate during the inventory period for the NFI plots and future climates for the populations and the species range), and their climatic differences between reference and inventory/future climate.

Code

##########
# DATASETS
##########

# We use the following datasets
  # clim_dfs which contains reference and future climates at the location of the populations
  # nfi_dfs which contains reference and survey climates at the location of the nfi plots
  # clim_ref_df which contains reference climates across the species distribution


# we also extract future climates across the species distribution (5 GCMs)
clim_fut_rast <- lapply(names(clim_dfs$clim_pred), function(gcm){ # for each GCM
  
# Rasters with future climates
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
clim_fut_df <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>% 
  mask(range_buffer) %>% 
  rasterToPoints() %>% # extract the climatic values at each spatial points
  as.data.frame()

# Scale the climatic variables with the initial scaling parameters
for(i in clim_var){clim_fut_df[,i] <- (clim_fut_df[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}

clim_fut_df
})
names(clim_fut_rast) <- names(clim_dfs$clim_pred)

Reference climates

Code

######################################
# Density plots of reference climates
######################################

lapply(clim_var, function(x){

   ## Legend title
 var_name <- extract_climatedt_metadata(var_clim = x) %>% 
    mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
    pull(var_legend)
 
  # Climatic values at the location of the populations
  tibble(group="Populations",value=clim_dfs$clim_ref[[x]]) %>% 
    bind_rows(tibble(group="NFI plots", value=nfi_dfs$clim_ref[[x]])) %>% 
    bind_rows(tibble(group="Species distribution", value=clim_ref_df$bio1)) %>% 
    
    ggplot(aes(x=value, group=group, fill=group)) +
    xlab(var_name) +
    geom_density(adjust=1.5, alpha=.4) +
    theme_bw() + 
    ggtitle("Reference climates (1901-1950)") +
    theme(legend.position = c(0.2,0.8),
          legend.title = element_blank())
})

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

Future climates

Code

lapply(clim_var, function(x){
  
  ## Legend title
  var_name <- extract_climatedt_metadata(var_clim = x) %>% 
    mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
    pull(var_legend)
 
  # Climatic values at the location of the populations
 
lapply(names(clim_dfs$clim_pred), function(gcm){
  
  tibble(group="Populations",
         value=clim_dfs$clim_pred[[gcm]][[x]]) %>% 
    bind_rows(tibble(group="NFI plots", 
                     value=nfi_dfs$clim_pred[[x]])) %>% 
    bind_rows(tibble(group="Species distribution", 
                     value=clim_fut_rast[[gcm]][[x]])) %>% 
    mutate(gcm=gcm)
  }) %>% 
  bind_rows() %>% 
  
  ggplot(aes(x=value, group=group, fill=group)) +
  xlab(var_name) +
  ylab("Density") +
  geom_density(adjust=1.5, alpha=.4) +
  theme_bw() +
  facet_wrap(~gcm)+
  ggtitle("Future climates and inventory climates") +
  theme(legend.position = c(0.85,0.2),
        legend.title = element_blank())

})

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

Climatic differences

We look at the climatic differences between reference climates and future/inventory climates.

Code

lapply(clim_var, function(x){
  
  ## Legend title
  var_name <- extract_climatedt_metadata(var_clim = x) %>% 
    mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
    pull(var_legend)

  lapply(names(clim_dfs$clim_pred), function(gcm){
  
  tibble(group="Populations",
         value_ref=clim_dfs$clim_ref[[x]],
         value_pred=clim_dfs$clim_pred[[gcm]][[x]]) %>% 
 bind_rows(tibble(group="NFI plots", 
        value_ref=nfi_dfs$clim_ref[[x]],
        value_pred=nfi_dfs$clim_pred[[x]])) %>% 
   bind_rows(tibble(group="Species distribution", 
        value_ref=clim_ref_df[[x]],
        value_pred=clim_fut_rast[[gcm]][[x]])) %>% 
   mutate(diff=value_pred-value_ref,
           gcm=gcm)
}) %>% bind_rows() %>% 
    
    ggplot(aes(x=diff, group=group, fill=group)) +
    xlab("Future/Inventory climate - reference climate") +
    ggtitle(paste0("Climatic differences - ",var_name)) +
    ylab("Density") +
    geom_density(adjust=1.5, alpha=.4) +
    theme_bw() +
    facet_wrap(~gcm) +
    theme(legend.position = c(0.85,0.2),
          legend.title = element_blank())
})

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

These last graphs are useful to understand what is going on.

bio4 is the most important variable to explain the gene-climate relationships and thus to predict GO. For bio4, the climatic differences at the location of the NFI plots have higher minimum and maximum values than climatic differences at the location of the studied populations or across the species range (under future climates). More precisely, some NFI plots experienced strong deviations in bio4 during the survey period, which are larger than the mean deviations expected under climate change.

This trend may seem surprising at first glance. However, it is important to remember that predictions of future climates are 20-year averages, whereas climates for NFI plots are calculated over shorter inventory periods. Consequently, NFI plot climates are likely to be subject to greater deviations from long-term climate averages, particularly for plots with shorter inventory periods.

Mapping

I generated maps to visualize the climatic differences but it does not help to visualize the higher climatic differences in some of the NFI plots for bio4.

Code

world <- ne_countries(scale = "medium", returnclass = "sf")

nfi_coord <- nfi_clim$clim_ref %>% dplyr::select(plotcode,longitude,latitude)

nfi_coord_spdf <- SpatialPointsDataFrame(coords = nfi_coord[, c("longitude", "latitude")],
                                         data = nfi_coord,
                                         proj4string = CRS("+proj=longlat +datum=WGS84"))


clim_nfi_fut <- lapply(names(clim_dfs$clim_pred), function(gcm){ # for each GCM
  
# Rasters with future climates
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))

df <- lapply(clim_var, function(x) paste0(path,"/",x,".tif")) %>% 
  raster::stack() %>% 
  extract(nfi_coord_spdf) %>% 
  as_tibble() %>% 
  bind_cols(nfi_coord) %>% 
  mutate(clim_period = paste0("clim_fut_",gcm))

}) %>% bind_rows

clim_nfi <- nfi_clim %>% 
  bind_rows(.id="clim_period") %>% 
  dplyr::select(clim_period, plotcode, longitude, latitude, any_of(clim_var)) %>% 
  bind_rows(clim_nfi_fut) %>% 
  pivot_longer(cols = any_of(clim_var), names_to = "var")

final_plots <- lapply(clim_var, function(x) {
  
   data <-   clim_nfi %>% dplyr::filter(var == x)
   
   max_value <- max(data$value, na.rm = T)
   min_value <- min(data$value, na.rm = T)

   
list_plots <-    lapply(unique(clim_nfi$clim_period), function(clim_period_i){
    
 subdata <-   data %>% dplyr::filter(clim_period == clim_period_i)
 
 
 ## Title of the graph
 if(clim_period_i == "clim_ref"){
   
   ggtitle <- "Reference climate (1901-1950)"
 } else if (clim_period_i == "clim_survey"){
   
   ggtitle <- "Climate during the inventory period"
 } else {
   ggtitle <- clim_period_i %>% str_sub(10,-1) %>% paste0("Future climate - ",.)
 }
 
 ## Legend title
 var_name <- extract_climatedt_metadata(var_clim = x) %>% 
    mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
    pull(var_legend)
 
 ## Plot
 ggplot() + 
   geom_sf(data = world, fill="gray98") + 
   theme_bw() +
   scale_x_continuous(limits = c(-10, 11)) +
   scale_y_continuous(limits = c(35, 51)) + 
   geom_point(data=subdata, aes(x=longitude,y=latitude,color=value), size=1, alpha = 0.8) + 
   xlab("") + ylab("") +
   ggtitle(ggtitle) +
   theme(legend.position = "right", # legend.position =c(0.85,0.2), legend.box.background = element_rect(fill = "white", linewidth=0.6, colour="gray")
          title = element_text(size=8),
          legend.text = element_text(size=10),
          legend.title = element_text(size=12)) +
      scale_colour_gradientn(colours=c("blue", "yellow", "red"), 
                              name = var_name,
                             na.value = "grey50",
                             limits=c(min_value,max_value))
  })
  
  
plot_legend <- get_legend(list_plots[[1]])

list_plots <- lapply(list_plots, function(p){
  p + theme(legend.position = "none")
})

list_plots$legend <- plot_legend
final_p <- plot_grid(plotlist = list_plots)

ggsave(final_p, filename=here(paste0("figs/GF/NFI_ClimaticMaps_",x,".pdf")), width=12, height=10)

final_p
})

lapply(final_plots, function(p) p)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

Relationship between GO and climatic differences

Finally, we look at the relationship between the climatic differences in the NFI plots (climate during the inventory period - climate during the reference period) and GO predictions.

We can see that the high GO values in the NFI plots are associated with strong deviations of bio4 from its values under the reference period. I think this is the reason why we obtained higher maximum GO values in the NFI plots than across the species range under future climates.

Code

lapply(clim_var, function(x){

   ## Legend title
 var_name <- extract_climatedt_metadata(var_clim = x) %>% 
    mutate(var_legend= paste0(description, " (", label,"; ",unit_symbol,")")) %>% 
    pull(var_legend)
 
tibble(group="NFI plots", 
       value_ref=nfi_dfs$clim_ref[[x]],
       value_pred=nfi_dfs$clim_pred[[x]],
       diff=value_pred - value_ref,
       go_pred=snp_sets$cand$go_nfi) %>% 
  ggplot(aes(x=diff, y=go_pred)) +
  geom_point(alpha=0.5) +
  xlab("Climate during the survey period - Reference climate") +
  ylab("Genomic offset predictions with GF and the candidate SNPs") +
  ggtitle(var_name) +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
  theme_bw()

})

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

Validation - Common gardens

Code

cg_clim <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>%  dplyr::select(cg,any_of(clim_var))
cg_coord <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,contains("ude"))
cg_names <- unique(cg_coord$cg)

# Generate scaled climatic datasets with climatic data at the location of the populations and at the location of the common gardens
cg_dfs <- generate_scaled_clim_datasets(clim_var, clim_pred = cg_clim)
  
# Predict genomic offset of each population when transplanted in the climate of the common gardens
snp_sets <- lapply(snp_sets, function(snp_set){

ref_pred <- predict(snp_set$gf_mod, as.data.frame(cg_dfs$clim_ref[,clim_var])) # predictions under reference-period climates
fut_pred <- predict(snp_set$gf_mod, as.data.frame(cg_dfs$clim_pred[,clim_var])) # predictions under climates during survey period

snp_set$go_cg <- lapply(1:nrow(ref_pred), function(x, ref_pred, fut_pred){
    as.numeric(pdist(ref_pred[x,],  fut_pred)@dist)}, fut_pred=fut_pred, ref_pred=ref_pred) %>% 
  setNames(cg_dfs$clim_ref[["pop"]]) %>% 
  as.data.frame() %>% 
  t() %>% 
  as.data.frame() %>% 
  set_colnames(cg_dfs$clim_pred[["cg"]]) %>% 
  rownames_to_column(var="pop") %>% 
  as_tibble()

return(snp_set)
})

# Map genomic offset predictions at the locations of the populations
go_maps_cg <- lapply(cg_names, function(cg_name){

# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(snp_set) snp_set$go_cg[[cg_name]]) %>%  unlist() %>% range()
go_limits[[1]] <- 0

p <- lapply(snp_sets, function(x){
            
 p <- make_go_map(dfcoord=pop_coord,
              snp_set = x,
              point_size = 3,
              type="CG",
              go_limits = go_limits,
              cg_name=cg_name,
              cg_coord=cg_coord)
 
 
 # We save the figures for the Supplementary Information
 if(x$set_code=="cand") {
   p_SI <- p + theme(plot.title = element_blank(),
                     legend.position = "none")
   
   p_SI %>% 
   ggsave(filename = here(paste0("figs/GF/GOmap_",x$set_code,"_",cg_name,"_SI.pdf")), device = "pdf",width=5,height=5)}
  
  if(x$set_code=="control") {
   p_SI <- p + theme(plot.title = element_blank())
   
   p_SI %>% 
   ggsave(filename = here(paste0("figs/GF/GOmap_",x$set_code,"_",cg_name,"_SI.pdf")), device = "pdf",width=7,height=5)}

 p
 
 })

plot_grid(p[[1]],p[[2]],p[[3]],nrow=1)
  
}) %>% setNames(cg_names)

pdf(here("figs/GF/GOmaps_CGs.pdf"), width=15,height=6)
lapply(go_maps_cg, function(x) x)
dev.off()

# show maps
lapply(go_maps_cg, function(x) x)

We look at the correlation across the different genomic offset predictions in the common gardens, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.

Code

par(mfrow=c(1,3))

lapply(cg_names, function(cg_name){

lapply(snp_sets, function(x) x$go_cg) %>% 

lapply(function(set){
  set[[cg_name]]
}) %>% 
    as_tibble() %>% 
    cor() %>% 
    corrplot(method = 'number',type = 'lower', 
             diag = FALSE,mar=c(0,0,2,0),
             title=str_to_title(cg_name),
             number.cex=2,tl.cex=1.5)

})

Let’s save the genomic offset predictions for comparison with the other methods.

Code

snp_sets %>% saveRDS(file=here("outputs/GF/go_predictions.rds"))

Session information

Code

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.3 (2024-02-29)
 os       Ubuntu 22.04.4 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2024-03-29
 pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package           * version date (UTC) lib source
 cachem              1.0.8   2023-05-01 [1] CRAN (R 4.3.2)
 cellranger          1.1.0   2016-07-27 [1] CRAN (R 4.3.2)
 class               7.3-22  2023-05-03 [4] CRAN (R 4.3.1)
 classInt            0.4-10  2023-09-05 [1] CRAN (R 4.3.2)
 cli                 3.6.2   2023-12-11 [1] CRAN (R 4.3.2)
 codetools           0.2-19  2023-02-01 [4] CRAN (R 4.2.2)
 colorspace          2.1-0   2023-01-23 [1] CRAN (R 4.3.2)
 corrplot          * 0.92    2021-11-18 [1] CRAN (R 4.3.2)
 cowplot           * 1.1.2   2023-12-15 [1] CRAN (R 4.3.2)
 DBI                 1.2.0   2023-12-21 [1] CRAN (R 4.3.2)
 devtools            2.4.5   2022-10-11 [1] CRAN (R 4.3.2)
 digest              0.6.34  2024-01-11 [1] CRAN (R 4.3.2)
 dplyr             * 1.1.4   2023-11-17 [1] CRAN (R 4.3.2)
 e1071               1.7-14  2023-12-06 [1] CRAN (R 4.3.2)
 ellipsis            0.3.2   2021-04-29 [1] CRAN (R 4.3.2)
 evaluate            0.23    2023-11-01 [1] CRAN (R 4.3.2)
 extendedForest      1.6.2   2023-12-12 [1] R-Forge (R 4.3.2)
 fansi               1.0.6   2023-12-08 [1] CRAN (R 4.3.2)
 farver              2.1.1   2022-07-06 [1] CRAN (R 4.3.2)
 fastmap             1.1.1   2023-02-24 [1] CRAN (R 4.3.2)
 forcats           * 1.0.0   2023-01-29 [1] CRAN (R 4.3.2)
 fs                  1.6.3   2023-07-20 [1] CRAN (R 4.3.2)
 generics            0.1.3   2022-07-05 [1] CRAN (R 4.3.2)
 ggplot2           * 3.4.4   2023-10-12 [1] CRAN (R 4.3.2)
 glue                1.6.2   2022-02-24 [1] CRAN (R 4.3.2)
 gradientForest    * 0.1-37  2023-08-19 [1] R-Forge (R 4.3.2)
 gtable              0.3.4   2023-08-21 [1] CRAN (R 4.3.2)
 here              * 1.0.1   2020-12-13 [1] CRAN (R 4.3.2)
 highr               0.10    2022-12-22 [1] CRAN (R 4.3.2)
 hms                 1.1.3   2023-03-21 [1] CRAN (R 4.3.2)
 htmltools           0.5.7   2023-11-03 [1] CRAN (R 4.3.2)
 htmlwidgets         1.6.4   2023-12-06 [1] CRAN (R 4.3.2)
 httpuv              1.6.13  2023-12-06 [1] CRAN (R 4.3.2)
 httr                1.4.7   2023-08-15 [1] CRAN (R 4.3.2)
 jsonlite            1.8.8   2023-12-04 [1] CRAN (R 4.3.2)
 kableExtra        * 1.3.4   2021-02-20 [1] CRAN (R 4.3.2)
 KernSmooth          2.23-22 2023-07-10 [4] CRAN (R 4.3.1)
 knitr             * 1.45    2023-10-30 [1] CRAN (R 4.3.2)
 labeling            0.4.3   2023-08-29 [1] CRAN (R 4.3.2)
 later               1.3.2   2023-12-06 [1] CRAN (R 4.3.2)
 latex2exp         * 0.9.6   2022-11-28 [1] CRAN (R 4.3.2)
 lattice             0.22-5  2023-10-24 [4] CRAN (R 4.3.1)
 lifecycle           1.0.4   2023-11-07 [1] CRAN (R 4.3.2)
 lubridate         * 1.9.3   2023-09-27 [1] CRAN (R 4.3.2)
 magrittr          * 2.0.3   2022-03-30 [1] CRAN (R 4.3.2)
 Matrix              1.6-3   2023-11-14 [4] CRAN (R 4.3.2)
 memoise             2.0.1   2021-11-26 [1] CRAN (R 4.3.2)
 mgcv                1.9-1   2023-12-21 [4] CRAN (R 4.3.2)
 mime                0.12    2021-09-28 [1] CRAN (R 4.3.2)
 miniUI              0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2)
 munsell             0.5.0   2018-06-12 [1] CRAN (R 4.3.2)
 nlme                3.1-163 2023-08-09 [4] CRAN (R 4.3.1)
 pdist             * 1.2.1   2022-05-02 [1] CRAN (R 4.3.2)
 pillar              1.9.0   2023-03-22 [1] CRAN (R 4.3.2)
 pkgbuild            1.4.3   2023-12-10 [1] CRAN (R 4.3.2)
 pkgconfig           2.0.3   2019-09-22 [1] CRAN (R 4.3.2)
 pkgload             1.3.3   2023-09-22 [1] CRAN (R 4.3.2)
 profvis             0.3.8   2023-05-02 [1] CRAN (R 4.3.2)
 promises            1.2.1   2023-08-10 [1] CRAN (R 4.3.2)
 proxy               0.4-27  2022-06-09 [1] CRAN (R 4.3.2)
 purrr             * 1.0.2   2023-08-10 [1] CRAN (R 4.3.2)
 R6                  2.5.1   2021-08-19 [1] CRAN (R 4.3.2)
 ragg                1.2.7   2023-12-11 [1] CRAN (R 4.3.2)
 raster            * 3.6-26  2023-10-14 [1] CRAN (R 4.3.2)
 RColorBrewer      * 1.1-3   2022-04-03 [1] CRAN (R 4.3.2)
 Rcpp                1.0.11  2023-07-06 [1] CRAN (R 4.3.2)
 readr             * 2.1.4   2023-02-10 [1] CRAN (R 4.3.2)
 readxl            * 1.4.3   2023-07-06 [1] CRAN (R 4.3.2)
 remotes             2.4.2.1 2023-07-18 [1] CRAN (R 4.3.2)
 rlang               1.1.3   2024-01-10 [1] CRAN (R 4.3.2)
 rmarkdown           2.25    2023-09-18 [1] CRAN (R 4.3.2)
 rnaturalearth     * 1.0.1   2023-12-15 [1] CRAN (R 4.3.2)
 rnaturalearthdata * 0.1.0   2017-02-21 [1] CRAN (R 4.3.2)
 rprojroot           2.0.4   2023-11-05 [1] CRAN (R 4.3.2)
 rstudioapi          0.15.0  2023-07-07 [1] CRAN (R 4.3.2)
 rvest               1.0.3   2022-08-19 [1] CRAN (R 4.3.2)
 scales              1.3.0   2023-11-28 [1] CRAN (R 4.3.2)
 sessioninfo         1.2.2   2021-12-06 [1] CRAN (R 4.3.2)
 sf                * 1.0-15  2023-12-18 [1] CRAN (R 4.3.2)
 shiny               1.8.0   2023-11-17 [1] CRAN (R 4.3.2)
 sp                * 2.1-2   2023-11-26 [1] CRAN (R 4.3.2)
 stringi             1.8.3   2023-12-11 [1] CRAN (R 4.3.2)
 stringr           * 1.5.1   2023-11-14 [1] CRAN (R 4.3.2)
 svglite             2.1.3   2023-12-08 [1] CRAN (R 4.3.2)
 systemfonts         1.0.5   2023-10-09 [1] CRAN (R 4.3.2)
 terra               1.7-65  2023-12-15 [1] CRAN (R 4.3.2)
 textshaping         0.3.7   2023-10-09 [1] CRAN (R 4.3.2)
 tibble            * 3.2.1   2023-03-20 [1] CRAN (R 4.3.2)
 tidyr             * 1.3.0   2023-01-24 [1] CRAN (R 4.3.2)
 tidyselect          1.2.0   2022-10-10 [1] CRAN (R 4.3.2)
 tidyverse         * 2.0.0   2023-02-22 [1] CRAN (R 4.3.2)
 timechange          0.2.0   2023-01-11 [1] CRAN (R 4.3.2)
 tzdb                0.4.0   2023-05-12 [1] CRAN (R 4.3.2)
 units               0.8-5   2023-11-28 [1] CRAN (R 4.3.2)
 urlchecker          1.0.1   2021-11-30 [1] CRAN (R 4.3.2)
 usethis             2.2.2   2023-07-06 [1] CRAN (R 4.3.2)
 utf8                1.2.4   2023-10-22 [1] CRAN (R 4.3.2)
 vctrs               0.6.5   2023-12-01 [1] CRAN (R 4.3.2)
 viridisLite         0.4.2   2023-05-02 [1] CRAN (R 4.3.2)
 webshot             0.5.5   2023-06-26 [1] CRAN (R 4.3.2)
 withr               3.0.0   2024-01-16 [1] CRAN (R 4.3.2)
 xfun                0.41    2023-11-01 [1] CRAN (R 4.3.2)
 xml2                1.3.6   2023-12-04 [1] CRAN (R 4.3.2)
 xtable              1.8-4   2019-04-21 [1] CRAN (R 4.3.2)
 yaml                2.3.8   2023-12-11 [1] CRAN (R 4.3.2)

 [1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.3
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Fitzpatrick, Matthew C, Vikram E Chhatre, Raju Y Soolanayakanahally, and Stephen R Keller. 2021. “Experimental Support for Genomic Prediction of Climate Maladaptation Using the Machine Learning Approach Gradient Forests.” Molecular Ecology Resources.